Diagnostic analysis of Aeolus return signal¶

Introduction¶

The SeaFlect project objective is to explore the information on sea surface wind contained in Aeolus surface return. This project considers the return signal from the so-called surface range bin. Before this return signal can be processed to arrive at a surface wind product, a diagnostic analysis of the return signal is performed to better understand the characteristics of this dataset.

The analysis considers oligotrophic waters as these are radiometrically stable in space and time. Also the storm track regions in the southern indian and atlantic oceans are considered as in these regions strong winds are common.

Before we look at the data a few words on the process envisaged to support the activities, to better understand the presented diagnostic analysis.

Dataset¶

The Aeolus dataset considered here is the so-called Mie useful signal (L1B data product). This signal is the Accumulated Charged Coupled Device (ACCD), which is an engineering unitless parameter. the ACCD-counts is provided for 24 layers. The interfaces of these layers are determined by the delay time of the return signal, i.e. the time it takes between sending and receiving the laser pulse. The nominal stratification of these layers is as follows, starting from the bottom : one pure ocean range bin (range bin 24), a composite range bin which contains an ocean layer, the sea surface and an atmospheric layer (range bin 23) and then 22 pure atmospheric layers. Here we consider the three lowest range bins. The 24 ACCD counts are complemented with a 25 parameter which contains the measured background signal. This data is also used for the diagnostic analysis. In summary the nominal numbers of the range bins considered here are in the next table

Number	Description
22	Lowest pure atmospheric layer
23	Surface layer
24	Top most pure ocean layer
25	Background signal

The thickness of the range bins are configurable parameters and they have been changed in-orbit. The number of layers is a fixed number.

The ACCD data comes in two two spatial resolution. There is one low spatial resolution dataset which represents the ACCD counts over an area of 87 km, and there is one high spatial resolution dataset, which represents the ACCD counts over an area of 3 km. For the presented analysis the course spatial resolution dataset has been considered.

Problem formulation¶

The problem at hand is to derive from the observations the reflectance (in physical units) of the combined surface interface and ocean layer. This because it is argued that this parameter is sensitive to the wind speed. Under low windspeed conditions this reflectance is small and dominated by the reflectance of the ocean layer. With increasing windspeed the white cap fraction increases. As white caps are bright, an increasing white cap fraction will result in an increase of the return signal.

The Aeolus L1b datastream provides the Mie-useful signal which is the so-called ACCD counts, for 24 range bins. This data is appended with the background signal. A total 25 elements are available. The data consistered here represent the Aeolus observations averaged over an 87 km track. This to increase the SNR.

The ACCD counts are engineering units and hence are unitless.

Furthermore the data of interest are contained in the surface range bin. Nominally this is bin number 23. The surface range bin is a composite layer, and consists of an atmospheric layer, the sea surface and an ocean layer. The thickness of these layers are not fixed, as instrument parameters which determine this are changing during in-orbit.

To extract the combined sea surface - ocean layer reflectance from the surface range bin, requires a correction of the atmospheric contribution. This contribution is estimated from the return signal from the lowest atmospheric layer, after a correction of the thickness of the atmospheric layers.

This implies that the thickness of the atmospheric sub-layer in the surface range bin, and of the lowest atmospheric layer are accurately knowns.

Furthermore it is assumed that the field of view is free of clouds or aerosol content, as we are interested in the return signal from the ocean-surface. Any contribution to the return signal from clouds-aerosol will only complicates the analysis.

Definitions¶

To follow the discussion presented here a number of variables are used.

Variable	Meaning
$R_{aso}$	Reflectance of the surface range bin
$R_{so}$	Combined reflectance of the sea-surface and ocean sub-layer
$R_{a}$	Reflectance of an atmospheric layer
$\Upsilon_{aso}$	The Accumulated Charged Coupled Device (ACCD) counts of the surface range bin
$\zeta$	The altitude of the upper boundary of the range gate [m]
$\Gamma$	The range of a range gate [m]

A a number superscript indicates the specific range bin, the nominal numbers are presented in the previous table.

Objective¶

The objective of the processing chain is to derive $R_{so}$ from $R_{aso}$. This requires a transformation from $\Upsilon_{aso}$ to the reflectance values $R_{aso}$, and further remove the atmospheric contribution $R_{a}^{23}$. This is done using $R_{a}^{22}$ after correcting this reflectance for layer difference.

The correction methodology itself is described below, but in general terms is based on a simple conceptual model which involves a standard approach to describe atmospheric multiple scattering. This conceptual model (CM), makes use of the the transmission of an atmospheric layer ${\cal T}$ which is a function of the number of molecules in a layer, so in general $${\cal T} = F\left( N_{molecules}\right) = F\left( P, T\right)$$. For the moment a fixed values for ${\cal t}$ is assumed. The conceptual model is described in more detail later.

Region of Interest¶

The analysis is performed for a number of regions over the global oceans. These regions are characterised by a chlorophyl content which is relatively stable over time. In particular we consider the South Pacific gyre which contains the “cleanest” water on earth and the Southern Storm track regions. The latter is included to have a large dynamical range of observed wind speeds.

Period of interest¶

The analysis considers the Aeolus observations for a 12 month period starting from september 2019.

Outline¶

What follows next is a walk through the various variables in the level 1b product stream for the region “E” (the south pacific gyre) for the entire period. Where relevant some detail plots are presented to highlight a certain feature.

Use of the notebook¶

To facilitate an easier flow of the results, most of the input cells are hidden. But they can be viewed by clicking on the icon if more details are needed.

Basic parameters¶

Thickness¶

The thickness of the range bins is important as it is used to estimate the atmospheric component in the surface range bin. There are two ways to estimate the thickness, namely from the altitude parameter of the interfaces between the range gates, and from the range parameter. The altitude parameter represents the height of the interfaces with respect to a reference geoid. the range parameter is a height of the range gate to the satellite.

For instance the thickness of range bin 22 is calculated either through

\[\Delta_\zeta^{22} = \zeta^{22} - \zeta^{23}\]

or through

\[\Delta_\Gamma^{22} = \Gamma^{23} - \Gamma^{22}\]

The next plots show the altitude and range parameters for the period of interest.

# the altitude of the interfaces
region="E"
graphic=[]
period="AnnualCycle"
dataset=AeolusProductObject.data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Height levels {} {}".format(period, region),\
         tools="box_zoom,pan,reset,save,wheel_zoom",toolbar_location="below")
#basic signal
p.circle(time, dataset["L1Bmie"]["alt"][:,range_bin["atmosphere"]], \
         color="red", legend_label="atmosphere")
p.circle(time, dataset["L1Bmie"]["alt"][:,range_bin["surface"]], \
         color="green", legend_label="surface")
p.circle(time, dataset["L1Bmie"]["alt"][:,range_bin["ocean"]],\
         color="blue", legend_label="ocean")
p.line(time, dataset["L1Bmie"]["altdem"][:], color="black", legend_label="dem")
p.yaxis.axis_label = 'Heigth [m]'
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300,toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"alt_time_serie_{}.png".format(period))
    export_png(plot, filename=figure_name)   

region="E"
graphic=[]
period="AnnualCycle"
dataset=AeolusProductObject.data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Ocean levels {} {}".format(period, region),\
         tools="box_zoom,pan,reset,save,wheel_zoom")
#basic signal
p.circle(time, dataset["L1Bmie"]["alt"][:,range_bin["ocean"]],\
         color="blue", legend_label="ocean")
p.line(time, dataset["L1Bmie"]["altdem"][:], color="black", legend_label="dem")
p.yaxis.axis_label = 'Heigth [m]'
p.y_range = Range1d(-150,50)
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300,toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"ocean_alt_time_serie_{}.png".format(period))
    export_png(plot, filename=figure_name)   

Discussion of the height signal¶

It appears that there are several peculiar values in the altitude parameter. The data at full dynamical range shows several isolated observations where the values of the altitude is well above any expected value. Outside of these isolated periods, the data looks reasonable.

ocean altitude¶

Although a detailed view is possible using the zoom function at the upper right of the figure as separate figure is presented with only the bottom interface of the surface range bin and the height of the surface according to the DEM. This zoom of the ocean altitude indicates a saw-tooth behaviour of the bottom boundary of the surface return, which is a common feature resulting from a index problem during interpolation.It is not clear if indeed the altitude is calculated using an interpolation routine, but if it is, then there is a hint of an index issue.

Second issue is the height of the ocean altitude with respect to the dem. Meaning during the beginning of the period of interest the bottom of the surface range bin appears to be very close to the surface as indicated by the DEM. The altitude is derived from the time difference between receiving and sending out the laser pulse. There is an uncertainty in the recording, which translates to an uncertainty of 100m in the altitude. This means that it is not clear that during the beginning of the period of interest the surface is located in range bin 22 (the nominal bin number) or in bin 24.

location of Surface¶

The location of the surface according to dem parameter is fluctuating by several meters. How this parameter is derived is unclear, It does not appear a critical parameters as the thickness of the atmospheric sub-layer in the surface range bin is of the order of 100 - 200 meters.

Range¶

The next figure shows the range variable for the three range bins considered here and this parameter shows a lot of variation over the time of interest. On first sight this appears normal as a satellite orbit is never fixed. What is remarkable is if a detailed of the range of a single overpass over the region of interest is analysed. Zoom is possible through the tools at the upper right of the figure. After zooming in it appears that the data is organised in the vertical, that is during a single overpass over the region of interest the range is changing by 6 km. And this patterns repeats itself the following day.

region="E"
graphic=[]
period="AnnualCycle"
dataset=AeolusProductObject.data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Range {} {}".format(period, region),\
         tools="box_zoom, pan,reset,save,wheel_zoom")
#basic signal
p.circle(time, dataset["L1Bmie"]["range"][:,range_bin["atmosphere"]], \
         color="red", legend_label="atmosphere")
p.yaxis.axis_label = 'Heigth [m]'
graphic.append(p)
p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Range {} {}".format(period, region),\
         tools="box_zoom, pan,reset,save,wheel_zoom")
#basic signal
p.circle(time, dataset["L1Bmie"]["range"][:,range_bin["surface"]], \
         color="green", legend_label="surface")
p.yaxis.axis_label = 'Heigth [m]'
graphic.append(p)

p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Range {} {}".format(period, region),\
         tools="box_zoom, pan,reset,save,wheel_zoom")
p.circle(time, dataset["L1Bmie"]["range"][:,range_bin["ocean"]],\
         color="blue", legend_label="ocean")
p.yaxis.axis_label = 'Heigth [m]'
graphic.append(p)

plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"range_time_serie_{}.png".format(period))
    export_png(plot, filename=figure_name)   

Thickness of range bin 21 (atmospheric range bin)¶

The introduction indicates that the thickness of a range bin can be calculated in two ways. It appears that these two methods do not generate the same value.

The next figure shows the value of $\Delta_\zeta^{22}$ and $\Delta_\Gamma^{22}$. These values are constant, over time, while at the same time a systematic difference between the two values can be observed. The value calculated according to the range being larger than the value according to the altitude. This is a relatively critical element as the thickness of the atmospheric range bin determines the amount of correction applied to the surface range bin. Hence which parameter is the correct one to use is important for the quality of the results.

region="E"
graphic=[]
period="AnnualCycle"
dataset=AeolusProductObject.data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
         title="Thickness atmospheric range bin {} {}".format(period, region),\
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom")
#basic signal
d1 = dataset["L1Bmie"]["range"][:,range_bin["surface"]] - dataset["L1Bmie"]["range"][:,range_bin["atmosphere"]]
d2 = dataset["L1Bmie"]["alt"][:,range_bin["atmosphere"]] - dataset["L1Bmie"]["alt"][:,range_bin["surface"]]
p.circle(time, d1, color="red", legend_label="range thickness")
p.circle(time, d2, color="green", legend_label="alt thickness")
p.yaxis.axis_label = 'Heigth [m]'
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"thickness_time_serie_{}.png".format(period))
    export_png(plot, filename=figure_name)   

region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]    
t=figure(x_axis_type="datetime", width=800, height=150, \
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom",\
           title="ACCD {} {}".format(period, region))
t.circle(time, dataset["L1Bmie"]["sgnl"][:,range_bin["top"]], \
         color="black", legend_label="top")
t.yaxis.axis_label = 'ACCD counts [-]'
t.legend.location = "top_left"
graphic.append(t)


s=figure(x_axis_type="datetime", width=800, height=150, \
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom",\
         title="ACCD {} {}".format(period, region))
s.circle(time, dataset["L1Bmie"]["sgnl"][:,range_bin["atmosphere"]], \
         color="red", legend_label="atmosphere")
s.yaxis.axis_label = 'ACCD counts [-]'
s.legend.location = "top_left"
graphic.append(s)

p22=figure(x_axis_type="datetime", width=800, height=150, \
           tools="box_zoom,pan,crosshair, reset,save,wheel_zoom",\
           title="ACCD {} {}".format(period, region))
p22.circle(time, dataset["L1Bmie"]["sgnl"][:,range_bin["surface"]], \
           color="green", legend_label="surface")
p22.yaxis.axis_label = 'ACCD counts [-]'
p22.legend.location = "top_left"
graphic.append(p22)


p23=figure(x_axis_type="datetime", width=800, height=150, \
           tools="box_zoom,pan,crosshair, reset,save,wheel_zoom",\
           title="ACCD {} {}".format(period, region))
p23.circle(time, dataset["L1Bmie"]["sgnl"][:,range_bin["ocean"]], \
           color="blue", legend_label="ocean")
p23.yaxis.axis_label = 'ACCD counts [-]'
p23.legend.location = "top_left"
graphic.append(p23)

plot = gridplot(graphic, ncols=1, width=800, height=300,\
                toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"basic_accd_signal_{}_{}.png".format(period, region))
    export_png(plot, filename=figure_name)   

Basic Signal¶

The above figures show the basic $\Upsilon$ as funtion of time for the 4 range bins. The three bins of interest and complemented with the top most range bin.

Return signal from upper most range bin¶

The return signal from the upper most range bin (in black) appears to have a dynamical range between 1000 and 3000 units. There are a few large excursions. There seems to be a trend to smaller return signals at the end of the time scale considered here. Although there is a link between the magnitude of the return signal and the thickness of the range bin, this thickness is relatively stable in time as we can see from previous plots, while there appears a linear trend. This is most likely linked to cross talk. As there is little mie return signal from the highest atmospheric layers, the return signal might be dominated by the rayleigh detector. It is known that the sensitivity is decreasing over time which could explain the trend.

Return signal from atmospheric and surface range bins¶

There are several interesting features before april 2020, but after this date, the return signals are relatively stable. No obvious trend is observed in either of the return signals. There is an indication that atmospheric signal is more varying towards end of the time domain, the ocean signal the opposite is observed.

Note¶

The reason the data is not normalised by thickness is that the surface range bin is a composite layer and hence such a normalisation does not automatically facilitate an easier comparison to other normalised return signals.

Background signal¶

The background signal shown below is the measured reflected solar irradiance. It is measured at a specific part of the detector array called the imaging zone. The solar backscattered signal is a function of the solar cycle, and hence it has an annual cyclus. It should be noted that the integration time of the backscatter (3750 $\mu$ s) is larger than the integration time of the lidar return signal (2.1 or 16.8 $\mu$s), hence the dynamical range is much larger than the return signals.

The backscatter has two curves, for the ascending and descending orbit.

The figure below the background signal shows the solar altitude as function of time. The solar altitude is calculated using a python routine based on the pvlib library. This library provides the azimuth and zenith angle, for a given location and time and height above the surface. For the calculations we assume the sea surface.

The solar altitude is calcated as 90 - zenith.

The azimuth angle is not presented as this provides no additional information.

If the solar altitude is smaller than a specific value, the sun is below the horizon and hence the amount of backscatter should be very small. As a threshold value of -7 or even a threshold of -16 can be used as the latter is defined as the astronomical twilight threshold.

During the periods April - September the sun is below the horizon, and the backscatter is relatively small.

What is remarkable from this figure is the sudden jump in solar altitude in spring and autumn. It is not clear if this is a bug in the diagnostic program (as the Aeolus time variable needs to be converted to a date-time string for the pvlib routine) or if the time variable itself has an jump.

region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]

#===========================================================
b=figure(x_axis_type="datetime", width=800, height=150, \
           title="Background {} {}".format(period, region),\
         tools="box_zoom, pan,reset,save,wheel_zoom")
b.circle(time, dataset["L1Bmie"]["sgnl"][:,range_bin["background"]], \
         color="brown", legend_label="background")
b.yaxis.axis_label = 'ACCD counts [-]'
b.legend.location = "top_left"
graphic.append(b)

s=figure(x_axis_type="datetime", plot_width=800, plot_height=150, \
           title="Solar Altitude {} {}".format(period, region),\
         tools="box_zoom, pan,reset,save,wheel_zoom")

s.y_range = Range1d(start=-50, end=50) 
s.extra_y_ranges = {"sunazi": Range1d(start=180, end=360) }

# s.circle(time, dataset["L1Bmie"]["sunazi"][:], color="green",\
#        y_range_name="sunazi")

s.circle(time, dataset["L1Bmie"]["sunalt"][:], color="blue")

x=[time[0], time[-1]]
y=[AeolusProductObject.sunalt_threshold, AeolusProductObject.sunalt_threshold]
s.line(x, y, color="black", line_width=6, \
       legend_label="Threshold")

y=[AeolusProductObject.astronomical_threshold, AeolusProductObject.astronomical_threshold]
s.line(x, y, color="purple", line_width=6, \
       legend_label="Astonomical Threshold")

#s.yaxis.axis_label = 'Altitude [Degree]'
## Adding the second axis to the plot. order is everything  
#s.add_layout(LinearAxis(y_range_name="sunazi", axis_label='Azimuth [Degree]'), 'right')

s.legend.location = "top_left"
graphic.append(s)

plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"background_solar_angles_{}_{}.png".format(period, region))
    export_png(plot, filename=figure_name)   

Variable	Meaning
\(R_{aso}\)	Reflectance of the surface range bin
\(R_{so}\)	Combined reflectance of the sea-surface and ocean sub-layer
\(R_{a}\)	Reflectance of an atmospheric layer
\(\Upsilon_{aso}\)	The Accumulated Charged Coupled Device (ACCD) counts of the surface range bin
\(\zeta\)	The altitude of the upper boundary of the range gate [m]
\(\Gamma\)	The range of a range gate [m]

SeaFlect

Diagnostic analysis of Aeolus return signal

Contents